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1.0  INTRODUCTION 


Until  recent  work  by  van  Dommelen  and  Shen  (1981,  1982a),  controversy 
existed  with  regard  to  the  possibility  that  a  singularity  occurs  at  a  finite 
time  and  position  in  an  unsteady  boundary  layer  evolving  under  prescribed 
pressure  distribution.  The  consequences  of  its  existence  are  large,  particu¬ 
larly  for  the  inviscid  flow.  It  seems  possible,  however,  that  previous 
uncertainties  have  arisen  from  inability  to  accurately  calculate  unsteady 
flows  with  large  reversal.  For  their  calculations,  van  Dommelen  and  Shen  made 
use  of  a  Lagrangian  formulation  of  the  two-dimensional  unsteady  boundary-layer 
equations  and  considered  a  circular  cylinder  started  impulsively  from  rest. 

They  confirmed  the  existence  of  a  singularity  and  the  accuracy  of  their  method. 
It  is  to  be  expected  that  the  singularity  will  also  exist  in  an  Eulerian 
formulation  but  many  previous  attempts,  for  example,  those  of  Belcher  et  al. 
(1971)  and  Telionis  and  Tsahalis  (1974)  have  failed  to  reveal  it.  Useful 
reviews  of  the  topic  have  been  provided  by  Riley  (1975),  Williams  (1977),  Shen 
(1978),  Telionis  (1979)  and  Cebeci  (1982). 

With  his  Eulerian  formulation,  Cebeci  (1979)  performed  calculations  for 
the  same  problem.  The  results  extended  to  x  =  2.8  at  which  time  they  were 
terminated  because  the  shear  layers  became  too  thick  to  be  accurately  repre¬ 
sented  by  the  numerical  scheme.  Also,  a  maximum  and  a  minimum  developed  in  the 
displacement  thickness  distribution  and  tended  to  cause  the  solutions  to 
oscillate.  Subsequent  calculations  with  an  improved  transformation  allowed 
solutions  up  to  x  =  3.0  but  the  nonmonotonic  variations  of  the  displacement 
thickness  remained,  Cebeci  (1982).  The  results  were  in  close  agreement  with 
those  of  van  Dommelen  and  Shen  (1982a)  up  to  around  t  =  2.75. 

Consideration  of  the  above  results  led  to  a  tentative  conclusion  that  the 
numerical  scheme  was  at  fault  and  that  previous  Eulerian  formulations  may  have 
had  similar  numerical  shortcomings.  Further  calculations  were  performed  to 
systematically  investigate  the  influence  of  step  lengths  in  time  and 
x-direction  and  confirmed  that  the  nonmonotonic  variation  of  displacement 
thickness  could  be  reduced  by  the  use  of  smaller  step  lengths  but  an  optimum 
relationship  between  At  and  Ax  could  not  be  determined  and  appeared  essential. 
These  calculations  were  performed  with  a  form  of  Keller's  Box  method  modified 
to  include  the  zig-zag  formulation  of  Krause  et  al.  (1968).  In  common  with 


the  often-used  Crank-Nicolson  method  (1947),  this  scheme  is  convenient,  partic¬ 
ularly  since  the  orientation  of  the  numerical  mesh  is  chosen  a  priori.  This 
advantage  has  a  corresponding  and  potentially  dangerous  drawback  in  the  pres¬ 
ence  of  large  reverse  flows.  The  precise  boundaries  of  the  solution  domain 
must  obey  the  Raetz  theory  (1957)  of  influence-regions  and  the  zones  of 
influence  and  dependence  may  not  be  properly  represented.  This  potential 
defect  can  be  overcome  only  by  ensuring  that  the  mesh  is  correctly  aligned 
with  the  flow  direction  and  thereby  allows  the  local  flow  properties  to  be 
correctly  influenced  by  its  neighbors. 

Here  we  describe  a  new  numerical  scheme  which  is  intelligent  in  the  sense 
that  it  examines  the  choice  of  grid  in  relation  to  the  magnitude  and  direction 
of  the  local  velocity  and  reaches  and  implements  a  decision  which  ensures  that 
Raetz 's  theory  is  obeyed.  It  is  applied  to  the  time  dependent  equations  of  the 
following  section  and  makes  use  of  a  form  of  the  characteristic  box  described 
in  Section  3.0.  The  results  for  an  impulsively  started  circular  cylinder  are 
presented  and  discussed  in  Section  4.0. 


2.0  BASIC  EQUATIONS 


We  consider  incompressible  unsteady  laminar  flow  over  a  circular  cylinder 
started  impulsively  from  rest.  The  governing  boundary-layer  equations  and 
their  boundary  conditions  for  this  flow  are  well  known,  see  for  example,  Cebeci 
(1979)  and  are  given  by 
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y  =  0,  u  =  v  =  0;  y  -*■  u  -*>  ug(x) 


(3) 


To  generate  the  initial  conditions  for  the  above  equations  at  t  =  0  and  to  put 
them  into  a  more  convenient  form  for  solution,  we  define  the  similarity  vari¬ 
able  n  and  a  dimensionless  stream  function  f(x,n,x)  by 
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where  ^  is  the  usual  definition  of  stream  function  which  satisfies  Eq.  (1), 
uQ  denotes  a  reference  velocity  2"ura,  Lisa  reference  length  equal  to  ra 
with  a  corresponding  to  the  cylinder  radius,  t  is  a  dimensionless  time  equal 
to  uQt/L,  ue  is  a  dimensionless  velocity,  ug/u0  and  x  =  x/L.  In  terms  of 
these  relations,  with  primes  denoting  differentiation  with  respect  to  n, 

Eqs.  (1)  to  (3)  can  be  written  as 


f "  1  +  -5-  f 11  +  T  — —  [1  -  ( f 1 ) 2  +  ff -]  =  x[y—  +  U  (f  '  —  -  f"  —  )]  (5) 

d  dx  31  e  3-  3X 

n  -  0 ,  f  =  f 1  =  0 ;  r>  =  rie,  f '  =  1  (6) 

These  variables  are  employed  only  in  the  interval  during  which  the  boundary 
layer  develops  rapidly.  At  higher  values  of  time,  the  exponential-1  ike  growth 
of  the  boundary-layer  thickness  with  time  near  the  rear  stagnation  point 
region,  is  represented  in  terms  of  the  dimensionless  variables  Y  and  F  defined 
by 


(7) 
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and  Eqs.  (l)-(3)  become 


0  =  (vu0L)1/,2eT_1ueFfx,Y,T) 


F"'  +  bnF"  +  b  —  [1  -  (F*  )2  +  FF"]  =  b  1^-  +  bu a[F'  -  F"  — ] 
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Here  b  =  e2^T”^  and  the  primes  denote  differentiation  with  respect  to  Y. 


Equations  (5)  and  (8)  can  also  be  expressed  in  the  common  form: 


f,M  +  b,nf "  +  b  —  [1  -  (f '  )2  +  ff-]  =  b~-  +  bu  (f '  -  f"  ^  ) 


where,  with  f  =  F  for  x  <  x* 


Ml  x  <_  x* 

b’"  t  >  t* 


T  1  T* 


X  >  Xj 


The  initial  conditions  at  t  =  0  can  be  obtained  from  Eq.  (10)  which  reduces  to 

f'"  +  j  nf"  =  0  (ID 

and  whose  solution,  subject  to  Eq.  (6),  is  given  by 

2 

f  =  n  erf  (£■)  +  [exp  (-  J-)  -  1]  (12) 


The  initial  conditions  along  the  (x,n)  plane  at  the  forward  and  rear  stagna¬ 
tion  points  can  also  be  obtained  from  Eq.  (10).  If  the  external  velocity 
distribution  is  represented  by 

u„  =  -  sin  (ttx  )  ( ‘ 

e  tt 

then  at  x  =  0  and  1,  Eq.  (10)  becomes 

f'M  +  b,n  f"  +  bX[l  -  (f')2  +  ff"]  =  b  C 

I  O  X 

where  X  =  1  for  x  =  0  and  X  =  -1  ror  x  =  1 . 
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3.0  NUMERICAL  METHOD 


The  solution  of  the  time-dependent  boundary-layer  equations  with  no  flow 
reversal  in  u  across  the  layer  is  routine  and  can  be  achieved  easily  with  the 
Keller  (1974)  or  Crank-Nicolson  method  (1947).  These  procedures  are  described 
in  several  references,  see  for  example,  Bradshaw  et  al.  (1981).  When  there  is 
flow  reversal,  however,  the  solution  of  the  time-dependent  boundary-layer 
equations  is  not  straightforward  and  requires  special  procedures  to  avoid  the 
numerical  instabilities.  The  zig-zag  formulation  of  Krause  et  al.  (1968)  helps 
to  remedy  the  difficulties  but  is  not  sufficient  in  all  cases,  as  we  shall 
discuss  later.  An  accurate  and  physically  sound  numerical  method  is  the 
characteristic  scheme  developed  for  three-dimensional  steady  boundary-layer 
flows  by  Cebeci  and  Stewartson  (1978)  which  is  based  on  the  solution  of 
governing  equations  along  local  streamlines.  It  allows  the  step  sizes  in  the 
t  and  x  directions  to  be  automatically  adjusted  to  ensure  that  the  region  of 
backflow  determined  by  the  local  streamlines  does  not  violate  Raetz's  regions 
of  influence.  A  brief  description  of  this  method  for  the  problem  of  the 
impulsively  started  cylinder  is  given  below. 


g'  =  w 


(17) 


and  write  Eqs.  (15)  and  (16)  and  their  boundary  conditions  as 
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dx  e  3x 
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w1  +  b,nw  +  0w  +  b  — —  (1  -  g^)  =  b  +  bu  g  — ^ 
1  d7  3t  e  - 


n  =  0, 


0  =  g  =  0; 


n  =  ne. 


1. 


(19) 

(20) 


The  solution  of  the  system  given  by  Eqs.  (17)  to  (20)  by  the  standard  or 
characteristic  box  methods  depends  on  how  the  difference  equations  are  written 
for  Eq.  (19);  the  remaining  equations  are  unchanged.  In  the  following  para¬ 
graphs  we  first  consider  Eqs.  (17)  and  (18)  and  show  how  the  difference  equa¬ 
tions  are  written. 

Consider  the  net  cube  shown  in  Fig.  1  and  denote  the  net  points  by 
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J 

The  difference  approximations  that  are  to  represent  Eq.  (17)  are  obtained 
by  averaging  about  the  midpoint  (x^*,  rn, 


.-1 i „i *  n  _ i , n \  , , 1  > " 

hj  (V  '  9j-l’  "  wM/2 


i  »n  ■ 


(22) 


where,  for  example. 


wl'i/2  ■  "2  +  «]:T>  ,23> 

The  difference  approximations  to  Eq.  (18)  are  obtained  by  centering  all 
quantities  except  0  at  the  center  of  the  cube  (x-j_i/2*  Tn- 1/2*  nj-l/2^  ^ 


6 


^  •srf/v  "-  JV'.n 


taking  the  values  of  each  variable,  say  q,  at  the  four  corners  of  the  box, 
that  is, 

qj-l/2*  =  2  (qj-l/2  +  qj-I/2)  =  J  (qj,n  +  ^j’1,n  +  Qjl"  +  q]l] *n)  (24a) 

However,  the  centering  of  0  is  done  by  writing  it  as 

In  terras  of  this  notation,  the  finite-difference  approximations  to  Eq. 

(18)  can  be  written  in  the  following  form: 


hTT  i~l/2,n-l/2 


hj  (9j  -9M}  '  (b 


9 j-1 /2 


=  (bu  )i“'(/2,n-l/2  -](  .  -  % 

e  i  yi-r 


where,  for  example, 

e  .  =  ei’V2,n-l/2^ 
J  J 


9j  =  7  (9j‘1/2,n  +  gj"172’""1). 


The  unknown  variables  in  Eq.  (25)  correspond  to  gVn  and  gi“l/2»n-V2  so 
that  when  a  solution  of  the  system  given  by  Eqs.  (17)  to  (20)  is  obtained,  g  is 
computed  at  (i,n,j)  and  6  at  (i-l/2,n-l/2, j).  This  modified  centering  pro¬ 
cedure  is  necessary  to  avoid  oscillations  due  to  the  use  of  the  continuity 
equation  in  the  form  given  by  Eq.  (15)  rather  than  the  use  of  the  stream  func¬ 
tion  which  allows  continuity  and  momentum  equations  to  be  expressed  as  a  third- 
order  equation. 

With  the  standard  box  scheme,  the  difference  approximations  corresponding 
to  Eq.  (19)  is  formulated  in  the  same  way  as  has  been  described  above.  With 
the  characteristic  box,  however,  an  alternative  procedure  is  followed  and  is 
described  below.  This  scheme  is  based  on  the  solution  of  the  governing  equa¬ 
tions  along  the  local  streamlines  defined  by 


If  we  denote  distance  in  this  direction  by  x,  and  the  angle  that  it  makes  with 
the  T-axis  by  8,  then  Eq.  (19)  can  be  written  as 


da. 


w 1  +  b,nw  +  0w  +  b  — -  (1  -  g2)  =  X 

1  dx  ^ 


(28) 


where 


X  =  b/1  +  (ueg)‘ 


(29a) 


1. 


8  =  tan  u 


(29b) 


The  finite-difference  approximations  to  Eq.  (28)  are  written  along  the 
streamline  direction  (see  Fig.  2).  The  resulting  expressions  are  lengthy  and 
to  illustrate  the  procedure,  we  consider  the  model  equation, 


w1  +  6w  =  X 


(30) 


With  the  notation  shown  in  Fig.  2,  the  difference  approximations  to  Eq.  (30)  at 
point  P  are 
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where  the  relation  between  9j_|/2  anc^  ^ose  values  of  6  centered  at  (i-1/2, 
n-1/2)  and  ( i -3/2 , n-1 /2 )  are 


i-3/2  i-1/2 
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(32) 


The  boundary  conditions  follow  from  Eqs.  (9)  and  can  be  written  as 


g  =  e  =  0; 
Jo  o  ’ 


9J  =  1 


(33) 


The  algebraic  system  given  by  Eqs.  (22),  (25)  and  (31)  together  with  the 
boundary  conditions  given  by  Eq.  (33)  is  nonlinear.  Linearization  is  achieved 


with  Newton's  method  and  the  equations  are  then  solved  by  the  block-elimination 
described,  for  example,  in  Bradshaw  et  al.  (1981). 

In  the  calculations  described  in  the  next  section,  the  standard  box  scheme 
has  been  used  for  all  situations  when  there  is  no  flow  reversal.  Where  a  calc¬ 
ulation  with  the  standard  box  reveals  a  flow  reversal  (u^  <  0)  further  itera¬ 
tions  at  that  location  make  use  of  the  characteristic  scheme  for  u.  <  0  and  the 

J 

standard  box  for  u-  >  0.  This  switch  from  one  scheme  to  another  continues  to 

Sj 

allow  quadratic  convergence  and  ensures  that,  provided  the  step  lengths  in  the 
t  and  x  directions  are  "properly"  selected,  the  numerical  instabilities  are 
avoided.  Further  details  are  provided  in  the  next  section. 


4.0  RESULTS 


t 

I 


It  is  convenient  to  present  the  results  which  we  have  obtained  using  the 
combination  of  standard  and  characteristic  box  schemes  and  then  to  discuss  them 
in  relation  to  previous  calculations.  Finally,  the  numerical  requirements  of 
the  present  results  are  discussed  together  with  the  implications  of  the  results 
for  future  calculations  of  more  realistic  unsteady  flows. 


Figures  3,4  and  5  display  the  variations  of  dimensionless  displacement 
thickness,  A*,  local  skin-friction  coefficient,  Cf  and  displacement 
velocity,  V  where  these  are  defined,  with  R.  =  u  L/v,  by 
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It  is  of  particular  note  that  the  displacement  thickness  is  close  to  monotonic 
with  the  small  maximum  and  minimum  for  t  =  3.1  at  which  the  calculations 
were  terminated.  The  previous  results  of  Cebeci  (1982)  are  also  shown  in  the 
figure  and  reveal  the  maxima  which  stemmed  from  the  use  of  a  numerical  scheme 
which  did  not  meet  the  requirements  imposed  by  the  Raetz  principle. 


The  distributions  of  local  skin-friction  coefficients  of  Fig.  4  show 
trends  which  are  similar  to  those  of  the  previous  results  but  with  differences 
in  magnitude  consistent  with  those  of  Fig.  3.  It  should  be  noted  that  the 
results  of  Figs.  3  and  4  are  identical  with  those  previously  obtained  up  to 
the  value  of  9  at  which  the  displacement  thickness  gradient  reaches  its 
maximum  and  for  values  of  t  less  than  around  2.75.  The  differences  for 
large  values  of  9  and  x  are  associated  with  the  numerical  procedure  and, 
in  particular,  with  its  ability  to  satisfy  the  Raetz  principle  as  is  discussed 
later. 


in 


The  dimensionless  displacement  velocity,  v  ,  is  shown  in  Fig.  5 
together  with  the  locus  of  points  corresponding  its  maxima  which  increases 
with  time  and  decreasing  angle.  At  t  =  3.0,  the  calculated  value  of  8  is 
111.5  and  corresponds  very  closely  to  that  determined  by  van  Dommelen  and  Shen 
who  terminated  their  calculations  at  this  time.  As  the  peak  in  the 
displacement  velocity  moves  upstream  with  increasing  time,  the  location  at 
which  the  skin-friction  coefficient  becomes  zero  also  moves  upstream  but  at  a 
slower  rate  and  towards  its  steady-state  value  of  105°,  Cebeci  and  Smith 
(1974).  It  may  be  conjectured  that  the  extrapolation  of  the  peak  displacement 
velocities  will  asymptote  to  this  steady-state  value  but  positive  confirmation 
will  require  excessive  computer  resources. 

Corresponding  velocity  profiles  are  presented  in  Fig.  6  for  t  =  2.75 
and  3.0.  The  general  trends  are  in  agreement  with  those  of  van  Dommelen  and 
Shen  and  that  the  quantitative  values  agree  closely  up  to  110°.  The 
discrepancies  at  larger  angles  are  probably  due  to  the  use  of  inappropriate 
time  steps  in  the  calculations.  The  previous  results  of  Cebeci  (1982)  tended 
to  agree  with  the  profiles  of  van  Dommelen  and  Shen  as  shown  in  figure  6a  and 
are  thought  to  be  inaccurate  as  is  discussed  in  the  following  section. 


5.0  CONCLUDING  REMARKS 


The  previous  calculations  of  Cebeci  (1982)  were  terminated  because  it 
proved  to  be  impossible  to  select  the  step  lengths  required  to  satisfy  the 
Raetz's  criteria  in  the  presence  of  large  flow  reversal.  Other  attempts  to 
make  use  of  Eulerian  formulations  to  calculate  unsteady  flows  with  large  back- 
flow  have  been  plagued  with  the  same  difficulty.  The  need  for  a  measure  which 
would  ensure  that  the  Raetz  principle  was  satisfied  has  been  met  in  the  present 
calculations.  Here  the  calculations  were  performed  at  a  given  time  and  with 
the  step  lengths  in  x  chosen  to  have  the  values  shown  in  Table  1.  As  can  be 
seen,  the  values  of  Ax  are  very  much  smaller  in  the  vicinity  of  singularity. 

The  solutions  were  iterated  at  each  x-station  until  a  convergence  criterion 

M 

based  on  the  wall  shear  parameter  fw  was  satisfied,  that  is, 

|(f;)V+I  -  tf“)v|  <  «,  (35) 
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where  6^  is  a  tolerance  parameter  which  was  set  equal  to  10  in  the 
calculations . 


The  keys  to  the  success  of  the  present  approach  lie  in  the  characteristic 
box  scheme  which  allows  the  orientation  of  the  finite-difference  mesh  to  vary 
across  the  shear  layer  and  the  procedure  for  the  automatic  selection  of  time 
steps  so  as  to  maintain  the  angle 


.  r. 

8  <  tan  p —  (36) 

Kn 

The  value  of  k  was  halved  until  this  condition  was  met.  The  resulting  val- 
n 

ues  of  kn  are  also  shown  in  Table  1  and  can  be  seen  to  become  extremely  small 
at  x  =  3.0.  The  present  calculations  which  made  use  of  increments  in  x,Y  and 
x  of  101,  161  and  435,  respectively,  could  have  been  extended  beyond  x  =  3.1 
but  at  considerable  expense  as  witnessed  by  the  small  and  decreasing  values  of 
kn*  The  values  of  kn  and  r^  shown  in  Table  1  were  subsequently  used  in  con¬ 
junction  with  the  zig-zag  scheme,  which  had  previously  failed,  Cebeci  (1982), 
to  permit  calculations  for  time  greater  than  x  >  2.75.  The  results  were 
found  to  be  identical  to  those  presented  here.  The  alternative  approach  of 
using  the  zig-zag  scheme  and  the  relationship  given  by  Eq.  (36)  was  not,  how¬ 
ever,  successful.  This  confirms  that  it  is  necessary  to  allow  the  flow  direc¬ 
tion  to  change  across  the  layer  by  conditions  determined  by  the  local  stream¬ 
lines  in  the  selection  of  k  . 


Table  1.  The  distribution  of  step  sizes  in  t  and  x. 


k  i 

r  • 

X 

n  | 

X 

T 

0  - 

1 

0.05 

1  o 

-  0.54 

0.02 

1  - 

1.5 

0.02 

1 

j  0.54 

+  0.57 

0.01 

1.5  - 

2.3 

0.01 

0.57 

-  0.58 

0.0025 

2.3  - 

2.73 

0.005 

0.58 

-  0.60 

0.0020 

2.73  - 

3.024 

0.002 

0.60 

-  0.612 

0.0015 

3.024  - 

3.1 

0.001 

0.612 

-  0.64 

0.0020 

0.64 

+  0.67 

0.0025 

0.67 

-  0.72 

0.01 

0.72 

-  1.0 

0.02 

It  is  clear  that  our  new  procedure  has  successfully  permitted  the  calcu¬ 
lation  of  the  flow  properties  for  the  unsteady  flow  associated  with  a  cylinder 
impulsively  started  from  rest.  The  large  reverse  flow  regions  found  with  this 
model  problem  occur  in  the  more  practical  application  of  oscillating  airfoils. 
Preliminary  work  has  been  carried  out  for  this  problem  and  has  confirmed  that 
this  is  so.  In  particular,  the  use  of  the  characteristic  box  scheme  together 
with  Eq.  (36)  led  to  solutions  which  approached  and  past  the  region  of  the 
singularity  without  numerical  difficulty  whereas  the  zig-zag  scheme  led  to 
solutions  which  oscillated  and  broke  down  in  the  same  manner  experienced  with 
the  cylinder. 

The  above  discussion  makes  it  clear  that  the  time  step  required  to  obtain 
results  at  the  larger  values  of  t  are  very  small  as  was  concluded  by  Ingham 
(1984).  For  this  reason,  the  present  calculations  were  terminated  at  x  =  3.1 
which  is  already  larger  than  that  of  previous  investigators.  Fi g-'  5  shows 

that  it  is  desirable  to  perform  calculations  at  higher  values  of  x  so  as  to 
confirm  the  conjecture  that  the  only  singularity  is  associated  with  the  steady- 
state  solution.  To  make  a  conclusive  judgment,  calculations  should  be  per¬ 
formed  up  to  x  =4.1  but,  as  Table  1  suggests,  the  required  time  steps  are 
likely  to  be  very  small.  The  time  required  to  obtain  results  in  the  range 
x  =  3.024  to  3.1,  which  corresponds  to  75  time-steps,  was  7  hours  on  a  CYBER 
175.  The  computer  time  likely  to  be  required  to  reach  x  =  4.1  is  clearly 
excessive. 
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Figure  3.  Variation  of  1 i-ens ion! ess  di spl acement  thicknes: 


Figure  5.  Variation  of  dimensionless  displacement  velocity.  The  insert  shows 
the  variation  of  t"e  location  of  maximum  displacement  velocity  ~m 
with  circles  indicating  the  computed  values,  the  dashed  lines  indicate 
the  linear  extrapolation  o*  m  and  the  solid  line  a  conjectured 
variation  of  to  steady  state. 
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